* Calibration sub-program for employment model (HW) with N-sectors (N) and
* variable hours (H) and intermediate inputs to production (I)

* Version "qfix_RP1" fixes recruiter productivity across industries at 37.5.


** Notes on calibration procedure

* Step 1: Load in raw data

* Step 2: Balance data to be consistent with SAM

* Step 3: Calibrate labor market variables and outer-nest of production function

* Step 4: Calibrate CES nest parameters

parameter
io0(ind,ind2)
iod0(ind,ind2)
iof0(ind,ind2)
c0(ind)
c_d0(ind)
c_f0(ind)
l0(ind)
lg0
ltot0
y0(ind)
e0(ind)
g0(ind)
g_d0(ind)
g_f0(ind)
yg0
exports0(ind)
;


$ONTEXT
$CALL GDXXRW.EXE pi22r.xlsx par=pi rng=A1:V2 cdim = 1
$CALL GDXXRW.EXE c22.xlsx par=c0 rng=A1:V2 cdim = 1
$CALL GDXXRW.EXE c_f22.xlsx par=c_f0 rng=A1:V2 cdim = 1
$CALL GDXXRW.EXE g22.xlsx par=g0 rng=A1:V2 cdim = 1
$CALL GDXXRW.EXE g_f22.xlsx par=g_f0 rng=A1:V2 cdim = 1
$CALL GDXXRW.EXE l22.xlsx par=l0 rng=A1:V2 cdim = 1
$CALL GDXXRW.EXE io22r.xlsx par=io0 rng=A1:W23 cdim = 1 rdim = 1
$CALL GDXXRW.EXE iof22.xlsx par=iof0 rng=A1:W23 cdim = 1 rdim = 1
$CALL GDXXRW.EXE e22.xlsx par=e0 rng=A1:V2 cdim = 1
$CALL GDXXRW.EXE ex22.xlsx par=exports0 rng=A1:V2 cdim = 1
$OFFTEXT


$GDXIN pi22r.gdx
$LOAD pi
$GDXIN

$GDXIN io22r.gdx
$LOAD io0
$GDXIN

$GDXIN iof22.gdx
$LOAD iof0
$GDXIN

$GDXIN c22.gdx
$LOAD c0
$GDXIN

$GDXIN c_f22.gdx
$LOAD c_f0
$GDXIN

$GDXIN l22.gdx
$LOAD l0
$GDXIN

$GDXIN g22.gdx
$LOAD g0
$GDXIN

$GDXIN g_f22.gdx
$LOAD g_f0
$GDXIN

$GDXIN e22.gdx
$LOAD e0
$GDXIN

$GDXIN ex22.gdx
$LOAD exports0
$GDXIN

c_d0(ind) = c0(ind) - c_f0(ind);
iod0(ind,ind2) = io0(ind,ind2) - iof0(ind,ind2);
lg0 = 1613.3;
ltot0 = sum(ind,l0(ind))+lg0;
g_d0(ind) = g0(ind) - g_f0(ind);

pig = pi('i22');
pi_f(ind) = pi(ind);
pig_f = pig;


y0(ind) = c_d0(ind) + g_d0(ind) + exports0(ind) + sum(ind2,iod0(ind,ind2));

** Scale consumption to get aggregate balance


parameter
ratiog
ratioc
;

ratiog = sum(ind,g0(ind))/sum(ind,c0(ind));

ratioc = sum(ind,l0(ind))/((1+ratiog)*sum(ind,c0(ind)));
c0(ind) = c0(ind)*ratioc;
c_d0(ind) = c_d0(ind)*ratioc;
c_f0(ind) = c_f0(ind)*ratioc;

g0(ind) = g0(ind)*ratioc;
g_f0(ind) = g_f0(ind)*ratioc;
g_d0(ind) = g0(ind) - g_f0(ind);

y0(ind) = c_d0(ind)+  g_d0(ind) + exports0(ind) + sum(ind2,iod0(ind,ind2));
yg0 = sum(ind,g0(ind))+lg0;

display
ratiog
ratioc
;

** Scale exports to get trade balance

parameter
imports0
ratiot
tbalance0
;

imports0(ind) = c_f0(ind) + g_f0(ind) + sum(ind2,iof0(ind,ind2));
ratiot = sum(ind,imports0(ind))/sum(ind,exports0(ind));
exports0(ind) = exports0(ind)*ratiot;
tbalance0 = sum(ind,exports0(ind)) - (sum(ind,c_f0(ind)) + sum(ind,g_f0(ind)) + sum(ind2,sum(ind,iof0(ind,ind2))));

display
ratiot
tbalance0
exports0
imports0
;

** Scale IO matrix to get industry balance

parameter
tfd
tfs
ex
extot
;

tfd(ind) = c_d0(ind)+  g_d0(ind) + exports0(ind) + sum(ind2,iod0(ind,ind2));
tfs(ind) = l0(ind) + sum(ind2,iod0(ind2,ind)+iof0(ind2,ind));
ex(ind) = tfd(ind)-tfs(ind);
extot = sum(ind,ex(ind));

display
tfd
tfs
ex
extot
;

$INCLUDE SCALE_TRADE.gms

y0(ind) = c_d0(ind) + g_d0(ind) + exports0(ind) + sum(ind2,iod0(ind,ind2));
tfd(ind) = y0(ind);
tfs(ind) = l0(ind) + sum(ind2,iod0(ind2,ind)+iof0(ind2,ind));
ex(ind) = tfd(ind)-tfs(ind);
io0(ind,ind2) = iod0(ind,ind2) + iof0(ind,ind2);

display
tfd
tfs
ex
;


* Calibrate emissions factor
parameter
f_y0 "oil/gas intensity"
etot0
e_y0
;


f_y0(secondary) = sum(oilgas,io0(oilgas,secondary))/y0(secondary);
*mu_e(ind) = e0(ind)/sum(ind2,(iod0(ind,ind2)+iof0(ind,ind2)));

mu_e(coal,elec) = 1.350/io0(coal,elec);
mu_e(coal,industrial) = .127/sum(industrial2,io0(coal,industrial2));
mu_e(coal,commercial) = .003/sum(commercial2,io0(coal,commercial2));
mu_e(oilgas,refine) = 2.296/(io0(oilgas,refine) - f_y0(refine)*exports0(refine) + f_y0(refine)*imports0(refine));
mu_e(oilgas,refine_n) = (1.473-.524)/(sum(refine_n2,io0(oilgas,refine_n2)) -sum(elec,io0(oilgas,elec)) - sum(ngd,f_y0(ngd)*exports0(ngd)) + sum(ngd,f_y0(ngd)*imports0(ngd)));
mu_e(oilgas,elec) = .524/(io0(oilgas,elec));

e0(ind) = sum(ie,mu_e(ie,ind)*io0(ie,ind));
etot0 = sum(ind,e0(ind)) - sum(oilgas,sum(secondary,mu_e(oilgas,secondary)*f_y0(secondary)*exports0(secondary))) + sum(oilgas,sum(secondary,mu_e(oilgas,secondary)*f_y0(secondary)*imports0(secondary)));

e0(ind) = e0(ind) + sum(secondary,sum(oilgas,mu_e(oilgas,secondary))*f_y0(secondary)*iof0(secondary,ind)) ;
e_y0(ind) = e0(ind)/y0(ind);

display
e_y0
f_y0
mu_e
etot0
;


* Scale data such that sum(c0) = 1;

parameter
ctot
;

ctot = sum(ind,c0(ind));
c0(ind) = c0(ind)/ctot;
c_d0(ind) = c_d0(ind)/ctot;
c_f0(ind) = c_f0(ind)/ctot;
y0(ind) = y0(ind)/ctot;
l0(ind) = l0(ind)/ctot;
lg0 = lg0/ctot;
g0(ind) = g0(ind)/ctot;
g_d0(ind) = g_d0(ind)/ctot;
g_f0(ind) = g_f0(ind)/ctot;
io0(ind,ind2) = io0(ind,ind2)/ctot;
iod0(ind,ind2) = iod0(ind,ind2)/ctot;
iof0(ind,ind2) = iof0(ind,ind2)/ctot;
e0(ind) = e0(ind)/ctot;
etot0 = etot0/ctot;
exports0(ind) = exports0(ind)/ctot;
imports0(ind) = imports0(ind)/ctot;
yg0 = yg0/ctot;

* Inpute Foreign Data

parameter
ratiof
sharec
shareg
shareio
l0_f
lg0_f
g0_f
g_d0_f
g_f0_f
io0_f
iod0_f
iof0_f
ctot_f
c0_f
c_d0_f
c_f0_f
imports0_f
exports0_f
y0_f
yg0_f
;

ratiof = 3;

sharec(ind) =  c0(ind)/(c0(ind)+g0(ind)+sum(ind2,io0(ind,ind2)));
shareg(ind) =  g0(ind)/(c0(ind)+g0(ind)+sum(ind2,io0(ind,ind2)));
shareio(ind,ind2) = io0(ind,ind2)/sum(ind3,io0(ind,ind3));


l0_f(ind) = ratiof*l0(ind);
lg0_f = ratiof*lg0;

g0_f(ind) = ratiof*g0(ind);
g_f0_f(ind) = exports0(ind)*shareg(ind);
g_d0_f(ind) = g0_f(ind) - g_f0_f(ind);

io0_f(ind,ind2) = ratiof*io0(ind,ind2);
iof0_f(ind,ind2) = exports0(ind)*(1-sharec(ind)-shareg(ind))*shareio(ind,ind2);
iod0_f(ind,ind2) = io0_f(ind,ind2) - iof0_f(ind,ind2);

ctot_f = ratiof*ctot;
c0_f(ind) = ratiof*c0(ind);
c_f0_f(ind) = exports0(ind)*sharec(ind);
c_d0_f(ind) =  c0_f(ind) - c_f0_f(ind);

imports0_f(ind) = sum(ind2,iof0_f(ind,ind2)) + c_f0_f(ind) + g_f0_f(ind);
exports0_f(ind) = sum(ind2,iof0(ind,ind2)) + c_f0(ind) + g_f0(ind);

y0_f(ind) = c_d0_f(ind) + g_d0_f(ind) + exports0_f(ind) + sum(ind2,iod0_f(ind,ind2));
yg0_f = sum(ind,g0_f(ind)) + lg0_f;

display
exports0
imports0_f
;

** Scale IO matrix to get industry balance

parameter
tfd_f
tfs_f
ex_f
extot_f
;

tfd_f(ind) = c_d0_f(ind)+  g_d0_f(ind) + exports0_f(ind) + sum(ind2,iod0_f(ind,ind2));
tfs_f(ind) = l0_f(ind) + sum(ind2,iod0_f(ind2,ind)+iof0_f(ind2,ind));
ex_f(ind) = tfd_f(ind)-tfs_f(ind);
extot_f = sum(ind,ex_f(ind));

display
tfd_f
tfs_f
ex_f
extot_f
;

$INCLUDE SCALE_TRADEF.gms

y0_f(ind) = c_d0_f(ind) + g_d0_f(ind) + exports0_f(ind) + sum(ind2,iod0_f(ind,ind2));
tfd_f(ind) = y0_f(ind);
tfs_f(ind) = l0_f(ind) + sum(ind2,iod0_f(ind2,ind)+iof0_f(ind2,ind));
ex_f(ind) = tfd_f(ind)-tfs_f(ind);
io0_f(ind,ind2) = iod0_f(ind,ind2) + iof0_f(ind,ind2);
yg0_f = sum(ind,g0_f(ind)) + lg0_f;

display
tfd_f
tfs_f
ex_f
;


** Calibrate labor market variables and production outer-nest

parameters
cbar_ss
cbar_f_ss
p_ss(ind)
p_f_ss
p_cbar_ss
p_cbar_f_ss
p_c_ss
p_c_f_ss
p_ex_ss
exch_ss
lambda_ss
lambda_f_ss
in_d(ind)
in_f(ind)
gtot0
gtot0_f
lambdag_ss
lambdag_f_ss
;




variables
ntot_ss
n_ss(ind)
ng_ss
u_ss

theta_i_ss
theta_g_ss
thetatot_ss
mu(ind)
mug
phi_i_ss
phi_g_ss


phitot_ss
v_n_ss(ind)
v_ng_ss


gammay0(ind)
alphap(ind)
q_ss(ind)
f_v_ss(ind)
f_l_ss(ind)
j_n_ss(ind)
vbar_ss(ind)
w_ss(ind)
h_ss(ind)


gammayg0
alphapg
qg_ss
f_vg_ss
f_lg_ss
cg_n_ss
vbarg_ss
wg_ss
hg_ss

psi
b0

ntot_f_ss
n_f_ss(ind)
ng_f_ss
theta_f_ss
thetag_f_ss
thetatot_f_ss
mu_f(ind)
mug_f
phi_f_ss
phig_f_ss


phitot_f_ss
v_n_f_ss(ind)
v_ng_f_ss

gammay0_f(ind)
alphap_f(ind)
q_f_ss(ind)
f_v_f_ss(ind)
f_l_f_ss(ind)
j_n_f_ss(ind)
vbar_f_ss(ind)
w_f_ss(ind)
h_f_ss(ind)

gammayg0_f
alphapg_f
qg_f_ss
f_vg_f_ss
f_lg_f_ss
cg_n_f_ss
vbarg_f_ss
wg_f_ss
hg_f_ss

psi_f
b0_f
;

equations
eqn_ntot_ss
eqn_u_ss

eqn_theta_i_ss
eqn_theta_g_ss
eqn_thetatot_ss
eqn_mu(ind)
eqn_mug
eqn_phi_i_ss
eqn_phi_g_ss

eqn_v_n_ss(ind)
eqn_v_ng_ss
eqn_phitot_ss

eqn_y_ss(ind)
eqn_firm_eom(ind)
eqn_foc_vbar_ss(ind)
eqn_foc_v_ss(ind)
eqn_foc_n_ss(ind)
eqn_f_l_ss(ind)
eqn_f_v_ss(ind)
eqn_nash_w_ss(ind)
eqn_nash_h_ss(ind)


eqn_yg_ss
eqn_gov_eom
eqn_foc_vbarg_ss
eqn_foc_vg_ss
eqn_foc_ng_ss
eqn_f_lg_ss
eqn_f_vg_ss
eqn_nash_wg_ss
eqn_nash_hg_ss

eqn_ntot_f_ss
eqn_theta_f_ss
eqn_thetag_f_ss
eqn_thetatot_f_ss
eqn_mu_f(ind)
eqn_mug_f
eqn_phi_f_ss
eqn_phig_f_ss

eqn_v_n_f_ss(ind)
eqn_v_ng_f_ss
eqn_phitot_f_ss

eqn_y_f_ss(ind)
eqn_firm_eom_f(ind)
eqn_foc_vbar_f_ss(ind)
eqn_foc_v_f_ss(ind)
eqn_foc_n_f_ss(ind)
eqn_f_l_f_ss(ind)
eqn_f_v_f_ss(ind)
eqn_nash_w_f_ss(ind)
eqn_nash_h_f_ss(ind)


eqn_yg_f_ss
eqn_gov_eom_f
eqn_foc_vbarg_f_ss
eqn_foc_vg_f_ss
eqn_foc_ng_f_ss
eqn_f_lg_f_ss
eqn_f_vg_f_ss
eqn_nash_wg_f_ss
eqn_nash_hg_f_ss
;


* Define fixed parameters

cbar_ss = sum(ind,c0(ind));
cbar_f_ss = sum(ind,c0_f(ind));

p_ss(ind) = 1;
p_cbar_ss = 1;
p_c_ss(ind) = 1;
lambdag_ss = 1;
lambdag_f_ss = 1;
exch_ss = 1;
p_f_ss(ind) = 1;
p_cbar_f_ss = 1;
p_c_f_ss(ind) = 1;

lambda_ss = (cbar_ss**(-sigma))/p_cbar_ss;
in_d(ind) = sum(ind2,io0(ind2,ind));
gtot0 = sum(ind,g0(ind));

lambda_f_ss = ((cbar_f_ss/ratiof)**(-sigma))/p_cbar_f_ss;
in_f(ind) = sum(ind2,io0_f(ind2,ind));
gtot0_f = sum(ind,g0_f(ind));


* Labor market equations

eqn_ntot_ss.. ntot_ss =e= sum(ind,n_ss(ind))+ng_ss;
eqn_u_ss.. u_ss =e= sum(ind2,pi(ind2)*n_ss(ind2))/(sum(ind,phi_i_ss(ind)) + phi_g_ss);

eqn_theta_i_ss(ind).. theta_i_ss(ind) =e= phi_i_ss(ind)/q_ss(ind);
eqn_theta_g_ss.. theta_g_ss =e= phi_g_ss/qg_ss;
eqn_thetatot_ss.. thetatot_ss =e= sum(ind,theta_i_ss(ind)) + theta_g_ss;
eqn_mu(ind).. mu(ind) =e= q_ss(ind)/(thetatot_ss**(-gamma));
eqn_mug.. mug =e= qg_ss/(thetatot_ss**(-gamma));
eqn_phi_i_ss(ind).. phi_i_ss(ind) =e= pi(ind)*(n_ss(ind)/(1-ntot_ss));
eqn_phi_g_ss.. phi_g_ss =e= pig*(ng_ss/(1-ntot_ss));



eqn_ntot_f_ss.. ntot_f_ss =e= sum(ind,n_f_ss(ind))+ng_f_ss;
eqn_theta_f_ss(ind).. theta_f_ss(ind) =e= phi_f_ss(ind)/q_f_ss(ind);
eqn_thetag_f_ss.. thetag_f_ss =e= phig_f_ss/qg_f_ss;
eqn_thetatot_f_ss.. thetatot_f_ss =e= (sum(ind,h_f_ss(ind)*vbar_f_ss(ind)*n_f_ss(ind))+hg_f_ss*vbarg_f_ss*ng_f_ss)/(ratiof-ntot_f_ss);
eqn_mu_f(ind).. mu_f(ind) =e=  q_f_ss(ind)*thetatot_f_ss**gamma;
eqn_mug_f.. mug_f =e= qg_f_ss*thetatot_f_ss**gamma;
eqn_phi_f_ss(ind).. phi_f_ss(ind) =e= pi_f(ind)*(n_f_ss(ind)/(ratiof-ntot_f_ss));
eqn_phig_f_ss.. phig_f_ss =e= pig_f*(ng_f_ss/(ratiof-ntot_f_ss));



* Household Equations

eqn_v_n_ss(ind).. v_n_ss(ind) =e= (lambda_ss*((1-tau_l0)*w_ss(ind)*h_ss(ind)-b0*p_cbar_ss) -psi*(h_ss(ind)**(1+chi)/(1+chi)) - phitot_ss)/(1-beta*(1-pi(ind)));
eqn_v_ng_ss.. v_ng_ss =e= (lambda_ss*((1-tau_l0)*wg_ss*hg_ss-b0*p_cbar_ss) -psi*(hg_ss**(1+chi)/(1+chi)) - phitot_ss)/(1-beta*(1-pig));
eqn_phitot_ss.. phitot_ss =e= beta*sum(ind,phi_i_ss(ind)*v_n_ss(ind))+beta*phi_g_ss*v_ng_ss;

eqn_v_n_f_ss(ind).. v_n_f_ss(ind) =e= (lambda_f_ss*((1-tau_l0)*w_f_ss(ind)*h_f_ss(ind)-b0_f*p_cbar_f_ss) -psi_f*(h_f_ss(ind)**(1+chi)/(1+chi)) - phitot_f_ss)/(1-beta*(1-pi_f(ind)));
eqn_v_ng_f_ss.. v_ng_f_ss =e= (lambda_f_ss*((1-tau_l0)*wg_f_ss*hg_f_ss-b0_f*p_cbar_f_ss) -psi_f*(hg_f_ss**(1+chi)/(1+chi)) - phitot_f_ss)/(1-beta*(1-pig_f));
eqn_phitot_f_ss.. phitot_f_ss =e= beta*sum(ind,phi_f_ss(ind)*v_n_f_ss(ind))+beta*phig_f_ss*v_ng_f_ss;


* Firm Equations

eqn_y_ss(ind)..  y0(ind) =e= gammay0(ind)*(alphap(ind)*(h_ss(ind)*n_ss(ind)*(1-vbar_ss(ind)))**rho(ind) + (1-alphap(ind))*in_d(ind)**rho(ind))**(1/rho(ind));
eqn_firm_eom(ind).. vbar_ss(ind) =e= pi(ind)/(q_ss(ind)*h_ss(ind));
eqn_foc_vbar_ss(ind).. q_ss(ind) =e= (p_ss(ind)*f_l_ss(ind))/(beta*j_n_ss(ind));
eqn_foc_v_ss(ind).. f_v_ss(ind) =e= 1;
eqn_foc_n_ss(ind).. f_l_ss(ind) =e= ((1-beta*(1-pi(ind)))*j_n_ss(ind) + (1+tau_p0)*w_ss(ind)*h_ss(ind))/(p_ss(ind)*h_ss(ind));
eqn_f_l_ss(ind).. f_l_ss(ind) =e= alphap(ind)*gammay0(ind)*( alphap(ind)* (h_ss(ind)*n_ss(ind)*(1-vbar_ss(ind)))**rho(ind) + (1-alphap(ind))*in_d(ind)**rho(ind) )**(1/rho(ind)-1)*(h_ss(ind)*n_ss(ind)*(1-vbar_ss(ind)))**(rho(ind)-1);
eqn_f_v_ss(ind).. f_v_ss(ind) =e= (1-alphap(ind))*gammay0(ind)*( alphap(ind)* (h_ss(ind)*n_ss(ind)*(1-vbar_ss(ind)))**rho(ind) + (1-alphap(ind))*in_d(ind)**rho(ind) )**(1/rho(ind)-1)*in_d(ind)**(rho(ind)-1);

eqn_y_f_ss(ind)..  y0_f(ind) =e= gammay0_f(ind)*(alphap_f(ind)*(h_f_ss(ind)*n_f_ss(ind)*(1-vbar_f_ss(ind)))**rho(ind) + (1-alphap_f(ind))*in_f(ind)**rho(ind))**(1/rho(ind));
eqn_firm_eom_f(ind).. vbar_f_ss(ind) =e= pi_f(ind)/(q_f_ss(ind)*h_f_ss(ind));
eqn_foc_vbar_f_ss(ind).. q_f_ss(ind) =e= (p_f_ss(ind)*f_l_f_ss(ind))/(beta*j_n_f_ss(ind));
eqn_foc_v_f_ss(ind).. f_v_f_ss(ind) =e= 1;
eqn_foc_n_f_ss(ind).. f_l_f_ss(ind) =e= ((1-beta*(1-pi_f(ind)))*j_n_f_ss(ind) + (1+tau_p0)*w_f_ss(ind)*h_f_ss(ind))/(p_f_ss(ind)*h_f_ss(ind));
eqn_f_l_f_ss(ind).. f_l_f_ss(ind) =e= alphap_f(ind)*gammay0_f(ind)*( alphap_f(ind)* (h_f_ss(ind)*n_f_ss(ind)*(1-vbar_f_ss(ind)))**rho(ind) + (1-alphap_f(ind))*in_f(ind)**rho(ind) )**(1/rho(ind)-1)*(h_f_ss(ind)*n_f_ss(ind)*(1-vbar_f_ss(ind)))**(rho(ind)-1);
eqn_f_v_f_ss(ind).. f_v_f_ss(ind) =e= (1-alphap_f(ind))*gammay0_f(ind)*( alphap_f(ind)* (h_f_ss(ind)*n_f_ss(ind)*(1-vbar_f_ss(ind)))**rho(ind) + (1-alphap_f(ind))*in_f(ind)**rho(ind) )**(1/rho(ind)-1)*in_f(ind)**(rho(ind)-1);


*Government Equations

eqn_yg_ss.. yg0 =e= gammayg0*(alphapg*(hg_ss*ng_ss*(1-vbarg_ss))**rhog + (1-alphapg)*gtot0**rhog)**(1/rhog);
eqn_gov_eom.. vbarg_ss =e= pig/(qg_ss*hg_ss);
eqn_foc_vbarg_ss.. qg_ss =e= -(lambdag_ss*f_lg_ss)/(beta*cg_n_ss);
eqn_foc_vg_ss.. f_vg_ss =e= 1/lambdag_ss;
eqn_foc_ng_ss.. f_lg_ss =e= ((beta*(1-pig)-1)*cg_n_ss + (1+tau_p0)*wg_ss*hg_ss)/(lambdag_ss*hg_ss);
eqn_f_lg_ss.. f_lg_ss =e= alphapg*gammayg0*( alphapg* (hg_ss*ng_ss*(1-vbarg_ss))**rhog + (1-alphapg)*gtot0**rhog )**(1/rhog-1)*(hg_ss*ng_ss*(1-vbarg_ss))**(rhog-1);
eqn_f_vg_ss.. f_vg_ss =e= (1-alphapg)*gammayg0*( alphapg* (hg_ss*ng_ss*(1-vbarg_ss))**rhog + (1-alphapg)*gtot0**rhog )**(1/rhog-1)*gtot0**(rhog-1);

eqn_yg_f_ss.. yg0_f =e= gammayg0_f*(alphapg_f*(hg_f_ss*ng_f_ss*(1-vbarg_f_ss))**rhog + (1-alphapg_f)*gtot0_f**rhog)**(1/rhog);
eqn_gov_eom_f.. vbarg_f_ss =e= pig_f/(qg_f_ss*hg_f_ss);
eqn_foc_vbarg_f_ss.. qg_f_ss =e= -(lambdag_f_ss*f_lg_f_ss)/(beta*cg_n_f_ss);
eqn_foc_vg_f_ss.. f_vg_f_ss =e= 1/lambdag_f_ss;
eqn_foc_ng_f_ss.. f_lg_f_ss =e= ((beta*(1-pig_f)-1)*cg_n_f_ss + (1+tau_p0)*wg_f_ss*hg_f_ss)/(lambdag_f_ss*hg_f_ss);
eqn_f_lg_f_ss.. f_lg_f_ss =e= alphapg_f*gammayg0_f*( alphapg_f* (hg_f_ss*ng_f_ss*(1-vbarg_f_ss))**rhog + (1-alphapg_f)*gtot0_f**rhog )**(1/rhog-1)*(hg_f_ss*ng_f_ss*(1-vbarg_f_ss))**(rhog-1);
eqn_f_vg_f_ss.. f_vg_f_ss =e= (1-alphapg_f)*gammayg0_f*( alphapg_f* (hg_f_ss*ng_f_ss*(1-vbarg_f_ss))**rhog + (1-alphapg_f)*gtot0_f**rhog )**(1/rhog-1)*gtot0_f**(rhog-1);

* Nash Bargaining Equations

eqn_nash_w_ss(ind).. j_n_ss(ind) =e= (tau/(1-tau))*((1+tau_p0)/(1-tau_l0))*(1/lambda_ss)*v_n_ss(ind);
eqn_nash_h_ss(ind).. h_ss(ind) =e= ((((1-tau_l0)/(1+tau_p0))*lambda_ss*(f_l_ss(ind)*p_ss(ind)))/psi)**(1/chi);

eqn_nash_wg_ss.. cg_n_ss =e= -(tau/(1-tau))*((1+tau_p0)/(1-tau_l0))*(1/lambda_ss)*v_ng_ss;
eqn_nash_hg_ss.. hg_ss =e= ((((1-tau_l0)/(1+tau_p0))*lambda_ss*(f_lg_ss*lambdag_ss))/psi)**(1/chi);

eqn_nash_w_f_ss(ind).. j_n_f_ss(ind) =e= (tau/(1-tau))*((1+tau_p0)/(1-tau_l0))*(1/lambda_f_ss)*v_n_f_ss(ind);
eqn_nash_h_f_ss(ind).. h_f_ss(ind) =e= ((((1-tau_l0)/(1+tau_p0))*lambda_f_ss*(f_l_f_ss(ind)*p_f_ss(ind)))/psi_f)**(1/chi);

eqn_nash_wg_f_ss.. cg_n_f_ss =e= -(tau/(1-tau))*((1+tau_p0)/(1-tau_l0))*(1/lambda_f_ss)*v_ng_f_ss;
eqn_nash_hg_f_ss.. hg_f_ss =e= ((((1-tau_l0)/(1+tau_p0))*lambda_f_ss*(f_lg_f_ss*lambdag_f_ss))/psi_f)**(1/chi);


* Initialization and fix variables

* Recruiter productivity fixed 

q_ss.fx(ind) = 37.5;
qg_ss.fx = 37.5;
q_f_ss.fx(ind) = 37.5;
qg_f_ss.fx = 37.5;

* Baseline unemployment set to 5 percent

ntot_ss.fx = 0.95;
ntot_f_ss.fx = 0.95*ratiof;

* Government hours fixed to 1/3 fraction of time

hg_ss.fx = 1/3;
hg_f_ss.fx = 1/3;

* Guess values of employment, disutility of work, unemployment benefits,
* and hours worked (in private industry)

v_n_ss.l(ind) = .106;
v_ng_ss.l = .106;
v_n_f_ss.l(ind) = .106;
v_ng_f_ss.l = .106;

psi.l = 7.9;
b0.l = .369;
psi_f.l = 7.9;
b0_f.l = .369;

h_ss.l(ind) = 1/3;
h_f_ss.l(ind) = 1/3;

* Use model equations to derive initialized values of remaining variables.

n_ss.l(ind) = l0(ind)*ntot_ss.l;
ng_ss.l = lg0*ntot_ss.l;
n_f_ss.l(ind) = n_ss.l(ind)*ratiof;
ng_f_ss.l = ng_ss.l*ratiof;

u_ss.l = (1-ntot_ss.l);



vbar_ss.l(ind) = pi(ind)/(q_ss.l(ind)*h_ss.l(ind));
vbarg_ss.l = pig/(qg_ss.l*hg_ss.l);

vbar_f_ss.l(ind) = pi_f(ind)/(q_f_ss.l(ind)*h_f_ss.l(ind));
vbarg_f_ss.l = pig_f/(qg_f_ss.l*hg_f_ss.l);



theta_i_ss.l(ind) = (vbar_ss.l(ind)*n_ss.l(ind)*h_ss.l(ind))/( u_ss.l);
theta_g_ss.l = hg_ss.l*vbarg_ss.l*ng_ss.l/( u_ss.l);
thetatot_ss.l = sum(ind,theta_i_ss.l(ind)) + theta_g_ss.l;
mu.l(ind) = q_ss.l(ind)/(thetatot_ss.l**(-gamma));
mug.l = qg_ss.l/(thetatot_ss.l**(-gamma) ) ;
phi_i_ss.l(ind) = mu.l(ind)*(theta_i_ss.l(ind)*thetatot_ss.l**(-gamma)  );
phi_g_ss.l = mug.l*(theta_g_ss.l*thetatot_ss.l**(-gamma));
phitot_ss.l = beta*sum(ind,phi_i_ss.l(ind)*v_n_ss.l(ind)) + beta*phi_g_ss.l*v_ng_ss.l;



phig_f_ss.l = pig_f*(ng_f_ss.l/(ratiof-ntot_f_ss.l));
phi_f_ss.l(ind) = pi_f(ind)*(n_f_ss.l(ind)/(ratiof-ntot_f_ss.l));
theta_f_ss.l(ind) = phi_f_ss.l(ind)/q_f_ss.l(ind);
thetag_f_ss.l = phig_f_ss.l/qg_f_ss.l;
thetatot_f_ss.l = (sum(ind,h_f_ss.l(ind)*vbar_f_ss.l(ind)*n_f_ss.l(ind))+hg_f_ss.l*vbarg_f_ss.l*ng_f_ss.l)/(ratiof-ntot_f_ss.l);
mu_f.l(ind) =  q_f_ss.l(ind)*thetatot_f_ss.l**gamma;
mug_f.l = qg_f_ss.l*thetatot_f_ss.l**gamma;
phitot_f_ss.l = beta*sum(ind,phi_f_ss.l(ind)*v_n_f_ss.l(ind))+beta*phig_f_ss.l*v_ng_f_ss.l;

w_ss.l(ind) =  (v_n_ss.l(ind)*(1-beta*(1-pi(ind))) + phitot_ss.l + psi.l*(h_ss.l(ind)**(1+chi))/(1+chi) + lambda_ss*p_cbar_ss*b0.l)/(lambda_ss*(1-tau_l0)*h_ss.l(ind));
wg_ss.l =  (v_ng_ss.l*(1-beta*(1-pig)) + phitot_ss.l + psi.l*(hg_ss.l**(1+chi))/(1+chi) + lambda_ss*p_cbar_ss*b0.l)/(lambda_ss*(1-tau_l0)*hg_ss.l);

w_f_ss.l(ind) =  (v_n_f_ss.l(ind)*(1-beta*(1-pi_f(ind))) + phitot_f_ss.l + psi_f.l*(h_f_ss.l(ind)**(1+chi))/(1+chi) + lambda_f_ss*p_cbar_f_ss*b0_f.l)/(lambda_f_ss*(1-tau_l0)*h_f_ss.l(ind));
wg_f_ss.l =  (v_ng_f_ss.l*(1-beta*(1-pig_f)) + phitot_f_ss.l + psi_f.l*(hg_f_ss.l**(1+chi))/(1+chi) + lambda_f_ss*p_cbar_f_ss*b0_f.l)/(lambda_f_ss*(1-tau_l0)*hg_f_ss.l);

j_n_ss.l(ind) = (tau/(1-tau))*((1+tau_p0)/(1-tau_l0))*(1/lambda_ss)*v_n_ss.l(ind);
cg_n_ss.l = -(tau/(1-tau))*((1+tau_p0)/(1-tau_l0))*(1/lambda_ss)*v_ng_ss.l;

j_n_f_ss.l(ind) = (tau/(1-tau))*((1+tau_p0)/(1-tau_l0))*(1/lambda_f_ss)*v_n_f_ss.l(ind);
cg_n_f_ss.l = -(tau/(1-tau))*((1+tau_p0)/(1-tau_l0))*(1/lambda_f_ss)*v_ng_f_ss.l;

f_l_ss.l(ind) = ((1-beta*(1-pi(ind)))*j_n_ss.l(ind) + (1+tau_p0)*w_ss.l(ind)*h_ss.l(ind))/(p_ss(ind)*h_ss.l(ind));
f_v_ss.l(ind) = 1;

f_lg_ss.l = ((beta*(1-pig)-1)*cg_n_ss.l + (1+tau_p0)*wg_ss.l*hg_ss.l)/(lambdag_ss*hg_ss.l);
f_vg_ss.l = 1;

f_l_f_ss.l(ind) = ((1-beta*(1-pi_f(ind)))*j_n_f_ss.l(ind) + (1+tau_p0)*w_f_ss.l(ind)*h_f_ss.l(ind))/(p_f_ss(ind)*h_f_ss.l(ind));
f_v_f_ss.l(ind) = 1;

f_lg_f_ss.l = ((beta*(1-pig_f)-1)*cg_n_f_ss.l + (1+tau_p0)*wg_f_ss.l*hg_f_ss.l)/(lambdag_f_ss*hg_f_ss.l);
f_vg_f_ss.l = 1;


* Note parameters x* and y* are used to derive initialized guesses for outer production function parameters

* (1/alpha) = x*/y* + 1 when f_v = 1 (f_l* = 1/x*)

parameter
x1
y1
x1g
y1g
x1_f
y1_f
x1g_f
y1g_f
;

x1(ind) = 1/(q_ss.l(ind)*beta*j_n_ss.l(ind));
y1(ind) = (in_d(ind)**(rho(ind)-1))/((h_ss.l(ind)*n_ss.l(ind)*(1-vbar_ss.l(ind)))**(rho(ind)-1));
x1g = -1/(qg_ss.l*beta*cg_n_ss.l);
y1g = (gtot0**(rhog-1))/((hg_ss.l*ng_ss.l*(1-vbarg_ss.l))**(rhog-1));
x1_f(ind) = 1/(q_f_ss.l(ind)*beta*j_n_f_ss.l(ind));
y1_f(ind) = (in_f(ind)**(rho(ind)-1))/((h_f_ss.l(ind)*n_f_ss.l(ind)*(1-vbar_f_ss.l(ind)))**(rho(ind)-1));
x1g_f = -1/(qg_f_ss.l*beta*cg_n_f_ss.l);
y1g_f = (gtot0_f**(rhog-1))/((hg_f_ss.l*ng_f_ss.l*(1-vbarg_f_ss.l))**(rhog-1));

alphap.l(ind) = y1(ind)/(x1(ind)+y1(ind));
gammay0.l(ind) = y0(ind)/(alphap.l(ind)*(h_ss.l(ind)*n_ss.l(ind)*(1-vbar_ss.l(ind)))**rho(ind) + (1-alphap.l(ind))*in_d(ind)**rho(ind))**(1/rho(ind)) ;
alphapg.l =  y1g/(x1g+y1g);
gammayg0.l = yg0/(alphapg.l*(hg_ss.l*ng_ss.l*(1-vbarg_ss.l))**rhog + (1-alphapg.l)*gtot0**rhog)**(1/rhog) ;
alphap_f.l(ind) = y1_f(ind)/(x1_f(ind)+y1_f(ind));
gammay0_f.l(ind) = y0_f(ind)/(alphap_f.l(ind)*(h_f_ss.l(ind)*n_f_ss.l(ind)*(1-vbar_f_ss.l(ind)))**rho(ind) + (1-alphap_f.l(ind))*in_f(ind)**rho(ind))**(1/rho(ind)) ;
alphapg_f.l =  y1g_f/(x1g_f+y1g_f);
gammayg0_f.l = yg0_f/(alphapg_f.l*(hg_f_ss.l*ng_f_ss.l*(1-vbarg_f_ss.l))**rhog + (1-alphapg_f.l)*gtot0_f**rhog)**(1/rhog) ;


* Define calibration model



model hwc_cal /
eqn_ntot_ss
eqn_u_ss

eqn_theta_i_ss
eqn_theta_g_ss
eqn_thetatot_ss
eqn_mu
eqn_mug
eqn_phi_i_ss
eqn_phi_g_ss

eqn_v_n_ss
eqn_v_ng_ss
eqn_phitot_ss

eqn_y_ss
eqn_firm_eom
eqn_foc_vbar_ss
eqn_foc_v_ss
eqn_foc_n_ss
eqn_f_l_ss
eqn_f_v_ss
eqn_nash_w_ss
eqn_nash_h_ss

eqn_yg_ss
eqn_gov_eom
eqn_foc_vbarg_ss
eqn_foc_vg_ss
eqn_foc_ng_ss
eqn_f_lg_ss
eqn_f_vg_ss
eqn_nash_wg_ss
eqn_nash_hg_ss

eqn_ntot_f_ss
eqn_theta_f_ss
eqn_thetag_f_ss
eqn_thetatot_f_ss
eqn_mu_f
eqn_mug_f
eqn_phi_f_ss
eqn_phig_f_ss

eqn_v_n_f_ss
eqn_v_ng_f_ss
eqn_phitot_f_ss

eqn_y_f_ss
eqn_firm_eom_f
eqn_foc_vbar_f_ss
eqn_foc_v_f_ss
eqn_foc_n_f_ss
eqn_f_l_f_ss
eqn_f_v_f_ss
eqn_nash_w_f_ss
eqn_nash_h_f_ss

eqn_yg_f_ss
eqn_gov_eom_f
eqn_foc_vbarg_f_ss
eqn_foc_vg_f_ss
eqn_foc_ng_f_ss
eqn_f_lg_f_ss
eqn_f_vg_f_ss
eqn_nash_wg_f_ss
eqn_nash_hg_f_ss

/
;

* Solve model

solve hwc_cal using mcp   ;

* Define post-solution parameters


parameter
utot_ss
utot_db
l_ss(ind)
lg_ss
tlump_ss
vlump_ss
grev_ss
m_ii
m_ig
m_gi
m_gg
m_from_i
m_from_g
m_ondiag

l_f_ss(ind)
lg_f_ss
tlump_f_ss
grev_f_ss

;


utot_ss = u_ss.l;
utot_db = 1-ntot_ss.l;
l_ss(ind) = h_ss.l(ind)*n_ss.l(ind);
lg_ss = hg_ss.l*ng_ss.l;
tlump_ss = (tau_l0+tau_p0)*(sum(ind,n_ss.l(ind)*w_ss.l(ind)*h_ss.l(ind)) + ng_ss.l*wg_ss.l*hg_ss.l) - (1-ntot_ss.l)*p_cbar_ss*b0.l - sum(ind,g0(ind)) - (1+tau_p0)*ng_ss.l*wg_ss.l*hg_ss.l;
vlump_ss = tlump_ss/(1-beta);
grev_ss = (tau_l0+tau_p0)*(sum(ind,n_ss.l(ind)*w_ss.l(ind)*h_ss.l(ind)) + ng_ss.l*wg_ss.l*hg_ss.l);




l_f_ss(ind) = h_f_ss.l(ind)*n_f_ss.l(ind);
lg_f_ss = hg_f_ss.l*ng_f_ss.l;
tlump_f_ss = (tau_l0+tau_p0)*(sum(ind,n_f_ss.l(ind)*w_f_ss.l(ind)*h_f_ss.l(ind)) + ng_f_ss.l*wg_f_ss.l*hg_f_ss.l) - (ratiof-ntot_f_ss.l)*p_cbar_f_ss*b0_f.l - sum(ind,g0_f(ind)) - (1+tau_p0)*ng_f_ss.l*wg_f_ss.l*hg_f_ss.l;
grev_f_ss = (tau_l0+tau_p0)*(sum(ind,n_f_ss.l(ind)*w_f_ss.l(ind)*h_f_ss.l(ind)) + ng_f_ss.l*wg_f_ss.l*hg_f_ss.l);

* Check for market clearance and balance

parameter
ex_d(ind)
prof_ss(ind)
bal_ss(ind)
balg_ss
bc_ss

prof_f_ss(ind)
bal_f_ss(ind)
balg_f_ss
bc_f_ss

ntot_db

;

* Market clearence and balance check
ex_d(ind) = y0(ind) - c_d0(ind) - g_d0(ind) - exports0(ind) - sum(ind2,iod0(ind,ind2));
prof_ss(ind) = p_ss(ind)*y0(ind) - in_d(ind) - (1+tau_p0)*w_ss.l(ind)*n_ss.l(ind)*h_ss.l(ind);
bal_ss(ind)  = (1+tau_p0)*w_ss.l(ind)*n_ss.l(ind)*h_ss.l(ind) + prof_ss(ind) - l0(ind);
balg_ss = (1+tau_p0)*ng_ss.l*wg_ss.l*hg_ss.l - lg0;
bc_ss = (1-tau_l0)*(sum(ind,n_ss.l(ind)*w_ss.l(ind)*h_ss.l(ind)) +ng_ss.l*wg_ss.l*hg_ss.l)  + (1-ntot_ss.l)*p_cbar_ss*b0.l + tlump_ss + sum(ind,prof_ss(ind)) - p_cbar_ss*cbar_ss;

ex_f(ind) = y0_f(ind) - c_d0_f(ind) - g_d0_f(ind) - imports0(ind) - sum(ind2,iod0_f(ind,ind2));
prof_f_ss(ind) = p_f_ss(ind)*y0_f(ind) - in_f(ind) - (1+tau_p0)*w_f_ss.l(ind)*n_f_ss.l(ind)*h_f_ss.l(ind);
bal_f_ss(ind)  = (1+tau_p0)*w_f_ss.l(ind)*n_f_ss.l(ind)*h_f_ss.l(ind) + prof_f_ss(ind) - l0_f(ind);
balg_f_ss = (1+tau_p0)*ng_f_ss.l*wg_f_ss.l*hg_f_ss.l - lg0_f;
bc_f_ss = (1-tau_l0)*(sum(ind,n_f_ss.l(ind)*w_f_ss.l(ind)*h_f_ss.l(ind)) +ng_f_ss.l*wg_f_ss.l*hg_f_ss.l)  + (ratiof-ntot_f_ss.l)*p_cbar_f_ss*b0_f.l + tlump_f_ss + sum(ind,prof_f_ss(ind)) - p_cbar_f_ss*cbar_f_ss;


ntot_db = sum(ind,n_ss.l(ind))+ng_ss.l ;


display
ex_d
prof_ss
bal_ss
balg_ss
bc_ss

ex_f
prof_f_ss
bal_f_ss
balg_f_ss
bc_f_ss


ntot_db
utot_db
utot_ss

;

* Display key calibration results
display
l0
v_n_ss.l
v_ng_ss.l
j_n_ss.l
cg_n_ss.l
q_ss.l
qg_ss.l
n_ss.l
ng_ss.l
ntot_ss.l
utot_ss
psi.l
b0.l
f_l_ss.l
gammay0.l
alphap.l
h_ss.l
hg_ss.l
w_ss.l
wg_ss.l
mu.l
mug.l
vbar_ss.l
vbarg_ss.l
theta_i_ss.l
theta_g_ss.l
thetatot_ss.l
prof_ss

l0_f
v_n_f_ss.l
v_ng_f_ss.l
j_n_f_ss.l
cg_n_f_ss.l
q_ss.l
qg_ss.l
n_f_ss.l
ng_f_ss.l
h_f_ss.l
hg_f_ss.l
w_f_ss.l
wg_f_ss.l
mu_f.l
mug_f.l
vbar_f_ss.l
vbarg_f_ss.l
theta_f_ss.l
thetag_f_ss.l
thetatot_f_ss.l
prof_f_ss

psi.l
psi_f.l
b0.l
b0_f.l

;


** Calibrate production inner nests

* US Domestic-foreign Composite

parameters
io_id0(ind,ind2)
io_if0(ind,ind2)
cost_df0(ind,ind2)
alphadf_d(ind,ind2)
alphadf_f(ind,ind2)
gammadf(ind,ind2)
;

loop(ind,
loop(ind2,
if(io0(ind,ind2) gt 0,

io_id0(ind,ind2) = iod0(ind,ind2)/io0(ind,ind2);
io_if0(ind,ind2) = iof0(ind,ind2)/io0(ind,ind2);
cost_df0(ind,ind2) = p_ss(ind)*io_id0(ind,ind2) + (p_f_ss(ind)/exch_ss)*io_if0(ind,ind2);


alphadf_d(ind,ind2) = (p_ss(ind)/cost_df0(ind,ind2))*(io_id0(ind,ind2))**(1/sig_df(ind));
alphadf_f(ind,ind2) = ((p_f_ss(ind)/exch_ss)/cost_df0(ind,ind2))*(io_if0(ind,ind2))**(1/sig_df(ind));

gammadf(ind,ind2) = alphadf_d(ind,ind2)+alphadf_f(ind,ind2);
alphadf_d(ind,ind2) = alphadf_d(ind,ind2)/gammadf(ind,ind2);
alphadf_f(ind,ind2) = alphadf_f(ind,ind2)/gammadf(ind,ind2);

else

io_id0(ind,ind2) = 1;
io_if0(ind,ind2) = 0;
cost_df0(ind,ind2) = 1;
alphadf_d(ind,ind2) = 1;
alphadf_f(ind,ind2) = 0;
gammadf(ind,ind2) = 1;
);
);
);

parameters
io_idg0(ind)
io_ifg0(ind)
cost_dfg0(ind)
alphadfg_d(ind)
alphadfg_f(ind)
gammadfg(ind)
;

loop(ind,
if(g0(ind) gt 0,

io_idg0(ind) = g_d0(ind)/g0(ind);
io_ifg0(ind) = g_f0(ind)/g0(ind);
cost_dfg0(ind) = p_ss(ind)*io_idg0(ind) + (p_f_ss(ind)/exch_ss)*io_ifg0(ind);


alphadfg_d(ind) = (p_ss(ind)/cost_dfg0(ind))*(io_idg0(ind))**(1/sig_df(ind));
alphadfg_f(ind) = ((p_f_ss(ind)/exch_ss)/cost_dfg0(ind))*(io_ifg0(ind))**(1/sig_df(ind));

gammadfg(ind) = alphadfg_d(ind)+alphadfg_f(ind);
alphadfg_d(ind) = alphadfg_d(ind)/gammadfg(ind);
alphadfg_f(ind) = alphadfg_f(ind)/gammadfg(ind);

else

io_idg0(ind) = 1;
io_ifg0(ind) = 0;
cost_dfg0(ind) = 1;
alphadfg_d(ind) = 1;
alphadfg_f(ind) = 0;
gammadfg(ind) = 1;
);
);


* ROW Domestic-foreign Composite

parameters
io_id0_f(ind,ind2)
io_if0_f(ind,ind2)
cost_df0_f(ind,ind2)
alphadf_d_f(ind,ind2)
alphadf_f_f(ind,ind2)
gammadf_f(ind,ind2)
;

loop(ind,
loop(ind2,
if(io0_f(ind,ind2) gt 0,

io_id0_f(ind,ind2) = iod0_f(ind,ind2)/io0_f(ind,ind2);
io_if0_f(ind,ind2) = iof0_f(ind,ind2)/io0_f(ind,ind2);
cost_df0_f(ind,ind2) = p_f_ss(ind)*io_id0_f(ind,ind2) + (p_ss(ind)*exch_ss)*io_if0_f(ind,ind2);


alphadf_d_f(ind,ind2) = (p_f_ss(ind)/cost_df0_f(ind,ind2))*(io_id0_f(ind,ind2))**(1/sig_df(ind));
alphadf_f_f(ind,ind2) = ((p_ss(ind)*exch_ss)/cost_df0_f(ind,ind2))*(io_if0_f(ind,ind2))**(1/sig_df(ind));

gammadf_f(ind,ind2) = alphadf_d_f(ind,ind2)+alphadf_f_f(ind,ind2);
alphadf_d_f(ind,ind2) = alphadf_d_f(ind,ind2)/gammadf_f(ind,ind2);
alphadf_f_f(ind,ind2) = alphadf_f_f(ind,ind2)/gammadf_f(ind,ind2);

else

io_id0_f(ind,ind2) = 1;
io_if0_f(ind,ind2) = 0;
cost_df0_f(ind,ind2) = 1;
alphadf_d_f(ind,ind2) = 1;
alphadf_f_f(ind,ind2) = 0;
gammadf_f(ind,ind2) = 1;
);
);
);

parameters
io_idg0_f(ind)
io_ifg0_f(ind)
cost_dfg0_f(ind)
alphadfg_d_f(ind)
alphadfg_f_f(ind)
gammadfg_f(ind)
;

loop(ind,
if(g0_f(ind) gt 0,

io_idg0_f(ind) = g_d0_f(ind)/g0_f(ind);
io_ifg0_f(ind) = g_f0_f(ind)/g0_f(ind);
cost_dfg0_f(ind) = p_f_ss(ind)*io_idg0_f(ind) + (p_ss(ind)*exch_ss)*io_ifg0_f(ind);


alphadfg_d_f(ind) = (p_f_ss(ind)/cost_dfg0_f(ind))*(io_idg0_f(ind))**(1/sig_df(ind));
alphadfg_f_f(ind) = ((p_ss(ind)*exch_ss)/cost_dfg0_f(ind))*(io_ifg0_f(ind))**(1/sig_df(ind));

gammadfg_f(ind) = alphadfg_d_f(ind)+alphadfg_f_f(ind);
alphadfg_d_f(ind) = alphadfg_d_f(ind)/gammadfg_f(ind);
alphadfg_f_f(ind) = alphadfg_f_f(ind)/gammadfg_f(ind);

else

io_idg0_f(ind) = 1;
io_ifg0_f(ind) = 0;
cost_dfg0_f(ind) = 1;
alphadfg_d_f(ind) = 1;
alphadfg_f_f(ind) = 0;
gammadfg_f(ind) = 1;
);
);

** Calibrate production inner nests

* Energy Composite

parameters
io_ie0(ie,ind)
cost_e0(ind)
alphae(ie,ind)
gammae(ind)
;

io_ie0(ie,ind) = io0(ie,ind)/sum(ie2,io0(ie2,ind));
cost_e0(ind) = sum(ie,cost_df0(ie,ind)*io_ie0(ie,ind));

alphae(ie,ind) = (cost_df0(ie,ind)/cost_e0(ind))*(io_ie0(ie,ind))**(1/sig_e(ind));
gammae(ind) = sum(ie,alphae(ie,ind));
alphae(ie,ind) = alphae(ie,ind)/gammae(ind);

parameters
io_ieg0(ie)
cost_eg0
alphaeg(ie)
gammaeg
;

io_ieg0(ie) = g0(ie)/sum(ie2,g0(ie2));
cost_eg0 = sum(ie,cost_dfg0(ie)*io_ieg0(ie));

alphaeg(ie) = (cost_dfg0(ie)/cost_eg0)*(io_ieg0(ie))**(1/sig_eg);
gammaeg = sum(ie,alphaeg(ie));
alphaeg(ie) = alphaeg(ie)/gammaeg;

parameters
io_ie0_f(ie,ind)
cost_e0_f(ind)
alphae_f(ie,ind)
gammae_f(ind)
;

io_ie0_f(ie,ind) = io0_f(ie,ind)/sum(ie2,io0_f(ie2,ind));
cost_e0_f(ind) = sum(ie,cost_df0_f(ie,ind)*io_ie0_f(ie,ind));

alphae_f(ie,ind) = (cost_df0_f(ie,ind)/cost_e0_f(ind))*(io_ie0_f(ie,ind))**(1/sig_e(ind));
gammae_f(ind) = sum(ie,alphae_f(ie,ind));
alphae_f(ie,ind) = alphae_f(ie,ind)/gammae_f(ind);

parameters
io_ieg0_f(ie)
cost_eg0_f
alphaeg_f(ie)
gammaeg_f
;

io_ieg0_f(ie) = g0_f(ie)/sum(ie2,g0_f(ie2));
cost_eg0_f = sum(ie,cost_dfg0_f(ie)*io_ieg0_f(ie));

alphaeg_f(ie) = (cost_dfg0_f(ie)/cost_eg0_f)*(io_ieg0_f(ie))**(1/sig_eg);
gammaeg_f = sum(ie,alphaeg_f(ie));
alphaeg_f(ie) = alphaeg_f(ie)/gammaeg_f;



* Material Composite

parameters
io_im0(im,ind)
cost_m0(ind)
alpham(im,ind)
gammam(ind)
;

io_im0(im,ind) = io0(im,ind)/sum(im2,io0(im2,ind));
cost_m0(ind) = sum(im,cost_df0(im,ind)*io_im0(im,ind));

alpham(im,ind) = (cost_df0(im,ind)/cost_m0(ind))*(io_im0(im,ind))**(1/sig_m(ind));
gammam(ind) = sum(im,alpham(im,ind));
alpham(im,ind) = alpham(im,ind)/gammam(ind);

parameters
io_img0(im)
cost_mg0
alphamg(im)
gammamg
;

io_img0(im) = g0(im)/sum(im2,g0(im2));
cost_mg0 = sum(im,cost_dfg0(im)*io_img0(im));

alphamg(im) = (cost_dfg0(im)/cost_mg0)*(io_img0(im))**(1/sig_mg);
gammamg = sum(im,alphamg(im));
alphamg(im) = alphamg(im)/gammamg;

parameters
io_im0_f(im,ind)
cost_m0_f(ind)
alpham_f(im,ind)
gammam_f(ind)
;

io_im0_f(im,ind) = io0_f(im,ind)/sum(im2,io0_f(im2,ind));
cost_m0_f(ind) = sum(im,cost_df0_f(im,ind)*io_im0_f(im,ind));

alpham_f(im,ind) = (cost_df0_f(im,ind)/cost_m0_f(ind))*(io_im0_f(im,ind))**(1/sig_m(ind));
gammam_f(ind) = sum(im,alpham_f(im,ind));
alpham_f(im,ind) = alpham_f(im,ind)/gammam_f(ind);

parameters
io_img0_f(im)
cost_mg0_f
alphamg_f(im)
gammamg_f
;

io_img0_f(im) = g0_f(im)/sum(im2,g0_f(im2));
cost_mg0_f = sum(im,cost_dfg0_f(im)*io_img0_f(im));

alphamg_f(im) = (cost_dfg0_f(im)/cost_mg0_f)*(io_img0_f(im))**(1/sig_mg);
gammamg_f = sum(im,alphamg_f(im));
alphamg_f(im) = alphamg_f(im)/gammamg_f;



* Intermediate Input Composite

parameters
vv0(ind)
io_ev0(ind)
io_mv0(ind)
cost_v0(ind)
alphav_e(ind)
alphav_m(ind)
gammav0(ind)
;


vv0(ind) = sum(ind2,io0(ind2,ind));
io_ev0(ind) = sum(ie,io0(ie,ind))/vv0(ind);
io_mv0(ind) = sum(im,io0(im,ind))/vv0(ind);
cost_v0(ind) = cost_e0(ind)*io_ev0(ind) + cost_m0(ind)*io_mv0(ind);

alphav_e(ind) = (cost_e0(ind)/cost_v0(ind))*(io_ev0(ind))**(1/sig_v(ind));
alphav_m(ind) = (cost_m0(ind)/cost_v0(ind))*(io_mv0(ind))**(1/sig_v(ind));

gammav0(ind) = alphav_e(ind)+alphav_m(ind);
alphav_e(ind) = alphav_e(ind)/gammav0(ind);
alphav_m(ind) = alphav_m(ind)/gammav0(ind);

parameters
vvg0
io_evg0
io_mvg0
cost_vg0
alphavg_e
alphavg_m
gammavg
;


vvg0 = sum(ind2,g0(ind2));
io_evg0 = sum(ie,g0(ie))/vvg0;
io_mvg0 = sum(im,g0(im))/vvg0;
cost_vg0 = cost_eg0*io_evg0 + cost_mg0*io_mvg0;

alphavg_e = (cost_eg0/cost_vg0)*(io_evg0)**(1/sig_vg);
alphavg_m = (cost_mg0/cost_vg0)*(io_mvg0)**(1/sig_vg);

gammavg = alphavg_e+alphavg_m;
alphavg_e = alphavg_e/gammavg;
alphavg_m = alphavg_m/gammavg;

parameters
vv0_f(ind)
io_ev0_f(ind)
io_mv0_f(ind)
cost_v0_f(ind)
alphav_e_f(ind)
alphav_m_f(ind)
gammav_f(ind)
;


vv0_f(ind) = sum(ind2,io0_f(ind2,ind));
io_ev0_f(ind) = sum(ie,io0_f(ie,ind))/vv0_f(ind);
io_mv0_f(ind) = sum(im,io0_f(im,ind))/vv0_f(ind);
cost_v0_f(ind) = cost_e0_f(ind)*io_ev0_f(ind) + cost_m0_f(ind)*io_mv0_f(ind);

alphav_e_f(ind) = (cost_e0_f(ind)/cost_v0_f(ind))*(io_ev0_f(ind))**(1/sig_v(ind));
alphav_m_f(ind) = (cost_m0_f(ind)/cost_v0_f(ind))*(io_mv0_f(ind))**(1/sig_v(ind));

gammav_f(ind) = alphav_e_f(ind)+alphav_m_f(ind);
alphav_e_f(ind) = alphav_e_f(ind)/gammav_f(ind);
alphav_m_f(ind) = alphav_m_f(ind)/gammav_f(ind);

parameters
vvg0_f
io_evg0_f
io_mvg0_f
cost_vg0_f
alphavg_e_f
alphavg_m_f
gammavg_f
;


vvg0_f = sum(ind2,g0_f(ind2));
io_evg0_f = sum(ie,g0_f(ie))/vvg0_f;
io_mvg0_f = sum(im,g0_f(im))/vvg0_f;
cost_vg0_f = cost_eg0_f*io_evg0_f + cost_mg0_f*io_mvg0_f;

alphavg_e_f = (cost_eg0_f/cost_vg0_f)*(io_evg0_f)**(1/sig_vg);
alphavg_m_f = (cost_mg0_f/cost_vg0_f)*(io_mvg0_f)**(1/sig_vg);

gammavg_f = alphavg_e_f+alphavg_m_f;
alphavg_e_f = alphavg_e_f/gammavg_f;
alphavg_m_f = alphavg_m_f/gammavg_f;



** Calibrate consumption aggregation function

* US Domestic-Foreign Consumption Good

parameter
alphac_d(ind)
alphac_f(ind)
gammacf(ind)
;



loop(ind,
if(c0(ind) gt 0,

alphac_d(ind) = (p_ss(ind)/p_c_ss(ind))*(c_d0(ind)/c0(ind))**(1/sig_c(ind));
alphac_f(ind) =  ((p_f_ss(ind)/exch_ss)/p_c_ss(ind))*(c_f0(ind)/c0(ind))**(1/sig_c(ind));

gammacf(ind) = alphac_d(ind)+alphac_f(ind);
alphac_d(ind) = alphac_d(ind)/gammacf(ind);
alphac_f(ind) = alphac_f(ind)/gammacf(ind);

else


alphac_d(ind) = 1;
alphac_f(ind) = 0;
gammacf(ind) = 1;
);
);

* ROW Domestic-Foreign Consumption Good

parameter
alphac_d_f(ind)
alphac_f_f(ind)
gammacf_f(ind)
;



loop(ind,
if(c0_f(ind) gt 0,

alphac_d_f(ind) = (p_f_ss(ind)/p_c_f_ss(ind))*(c_d0_f(ind)/c0_f(ind))**(1/sig_c_f(ind));
alphac_f_f(ind) =  ((p_ss(ind)*exch_ss)/p_c_f_ss(ind))*(c_f0_f(ind)/c0_f(ind))**(1/sig_c_f(ind));

gammacf_f(ind) = alphac_d_f(ind)+alphac_f_f(ind);
alphac_d_f(ind) = alphac_d_f(ind)/gammacf_f(ind);
alphac_f_f(ind) = alphac_f_f(ind)/gammacf_f(ind);

else


alphac_d_f(ind) = 1;
alphac_f_f(ind) = 0;
gammacf_f(ind) = 1;
);
);

* Consumption Bundle

parameter
alphac(ind)
gammac
;

alphac(ind) = (p_c_ss(ind)/p_cbar_ss)*(c0(ind)/cbar_ss)**(1/sigmac);
gammac = sum(ind,alphac(ind));
alphac(ind) = alphac(ind)/gammac;


parameter
alphac_ff(ind)
gammac_ff
;

alphac_ff(ind) = (p_c_f_ss(ind)/p_cbar_f_ss)*(c0_f(ind)/cbar_f_ss)**(1/sigmac);
gammac_ff = sum(ind,alphac_ff(ind));
alphac_ff(ind) = alphac_ff(ind)/gammac_ff;

* Used in the sensitivity analysis when unemployment benefit is proporational to after tax real wage;
parameter
b_rr
;

b_rr = (b0.l)/(( (sum(ind,(1-tau_l0)*w_ss.l(ind)*h_ss.l(ind)*n_ss.l(ind)) + (1-tau_l0)*wg_ss.l*hg_ss.l*ng_ss.l)/(sum(ind,n_ss.l(ind)) + ng_ss.l)  )/p_cbar_ss);

display
b_rr
;
